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ABSTRACT 


A.  A  computer  model  is  given  for  the  numerical  solution  of  the 
partial  differential  equations  governing  the  transport  and 
condensation  processes  of  clouds.  The  models  considered 
include  two  dimensional  warm  fog  seeding,  one  dimensional 
cumulus  cloud  formation,  and  three  dimensional  warm  fog 
modification  by  external  heat  sources.  The  finite  difference 
techniques  used  and  the  stability  criteria  are  discussed. 

The  method  of  calculating,  the  condensation/evaporation 
rates,  and  the  diffusion  parameters  are  also  given. 

B.  A  series  of  computer  programs,  data  reduction  techniques, 
and  data  analysis  procedures  are  developed  to  recover  per¬ 
tinent  data  from  various  recording  media  generated  on  AFCRL 
C-130  aircraft.  These  programs  and  procedures  are  developed 
for  the  PDP8/I  resident  at  AFCRL.  The  meteorological  data 
collected  is  reduced  into  a  usable  form  such  that  computer 
studies  relative  to  the  growth,  structure,  and  modification 
techniques  of  convective  clouds  can  be  made. 
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1.0  Introduction 


The  problems  that  were  considered  during  the  extent 
of  this  contract  were  to  develop  stable  numerical  techniques 
to  solve  the  non-linear  partial  differential  equations 
governing  the  .behavior  of  warm  fog  and  cumulus  clouds. 

These  techniques  must  be  valid  for  several  forms  of  the 
same  general  equations.  These  different  forms  would  be 
characterized  by  the  simplifying  assumptions  that  were  made 
to  reduce  the  general  equations  used  for  a  particular  model. 
Much  of  the  development  involved  numerical  experiments  to 
determine  the  sensitivity  of  the  models  to  the  numerical 
procedures.  X  major  portion  of  the  development  time  was 
spent  in  conference  with  the  Air  Force  to  insure  that  the 
physical  processes  were  accurately  represented  by  the  models. 

This  final  report  will  deal  with  the  most  successful 
experiments  only.  These  experiments  often  involved  the 
generation  of  computer  programs  to  validate  the  numerical 
procedures  and  the  simplifying  assumptions  used  in  applying 
them.  The  contract  considered  three  major  models,  the 
equations  for  which  were  provided  by  the  contract  monitor. 

1)  A  two  dimensional  Euleriar.  model  of  the  modifica¬ 
tion  of  warm  fog  by  hygroscopic  particle  seeding. 

2)  A  one  dimensional  warm  cumulus  model  involving 
detailed  calculations  of  the  condensation  and 
stochastic  coalescence  processes. 

3)  A  three  dimensional  Eulerian  model  of  warm  fog 
modification  by  the  addition  of  heat  at  the 
ground  level  from  mult:ple  point  and  line  sources. 

The  verification  of  the  models  considered  above  can  be 
made  by  collecting  data  during  the  actual  experiments.  Section 
3  deals  with  the  techniques  cf  reducing  the  data  on  various 


recording  media  collected  by  the  AFCRL  C-130  aircraft. 
The  techniques  of  data  reduction  are  presented  in  the 
form  of  procedures  that  also  represent  computer  programs 
for  the  POP-81  resident  of  AFCRL. 
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2.0  FOG  &  CLOUD  MODEL  SIMULATION 


2 . 1  Equations 


The  general  equations  used  to  approximate  the  warm  fog  and 
cloud  models  are  represented  by  the  diffusive  quantities  of 
momentum,  temperature,  vapor  and  particles.  The  equations 
of  motion  are  given  as: 


Ut  = 

Vt  =  -uvx 
wt  *  - “W X 


-  "uz  +  '  ?  Px 

-  »v2  ♦  (Vy’y  -  l  P} 

-  ww2  *  (Kmw2)2  +  B 


where  u,  v,  and  w  represent  the  velocities  In  the  r. ,  y  and 

r  directions  respecti vely ;  K  the  turbulent  diffusion  co¬ 
in 

efficient  for  momentum, B  the  buoyancy  term,,  and  P  the  pres¬ 
sure. 


The  temperature  equation  is: 

•t  ■  ■«',  -  v6y  ■  «,  *  +  + 

(V,1,  +  <KHez>z 


where  a  is  the  potential  temperature,  G  is  the  condensed 
water,  L  is  the  latent  heat,  Cp  is  the  specific  heat  and 
is  the  heat  diffusion  coefficient. 

The  mixing  ratio  equation  is: 

Mt  -  -UHx  -  v«y  -  wM2  -  f  ♦  (KeHx)x  ♦ 

(KEMyJy  +  (KEHz} z 

where  M  is  the  .mixing  ratio;  is  the  vapor  diffusion 
coeff  ic??'.nt. 


O',  the  modified  coefficient  of  molecular  diffusion  is 
given  as 


D’  =  ( . 001 5T  -  -18374)r  /  (r  +  .000486622) 


The  latent  heat,  L,  is  given  as 


L  =  750.2636  -  .56T 


(17) 


(18) 


k,  the  coefficient  of  thermal  conductivity  is  given  as 


k  =  (1.1 5628x1 0" 5  +  .\7xlO"6T)r/(r  +  1.929xl0~5) 


(19) 


The  gas  constant  for  water  vapor  is 


Rv  =  4.615x10 


(20) 


and  J  is  a  conversion  constant  4.18684x10  . 


The  diffusion  coefficients  Km,KH,  Kg,  and  Kp  are  evaluated 
as  functions  of  the  Richardson  number,  a  quantity  that 
measures  the  degree  of  air  turbulence.  The  Richardson 
number,  R..,  is  evaluated  as  follows: 


Rf  =  9.8e  /  (euz) 


(21) 


.!■ 


Ri  *  -2* 


fi^.25 

R.<-2. 


•2 .  ^  R,.  .  25 


(22) 


. . UiHaiMMatiuuil . . 


Two  quantities  used  to  directly  calculate  the  diffusion 
coefficients  are  $  and  a  which  are  functions  of  R-.  For 
the  unstable  condition  (R.<  0), 


$  =  (1  -  !5R.)"1/4 
a  =  1.35(1  -  9Ri)1,<4* 


(23) 


For  the  stable  condition*  R.  >  0, 


=  1  +  4. 7c 

a  «  *  /  (.74  +  4.7?) 


(25) 

(26) 


where 


.74  -  9.4R- 


(27) 


Having  $  and  a,  the  diffusion  coefficients  are 


Km  =  .  1 225’2  |u  ZH£\ 

KH  =  K  =  Kn  =  oK 
Hep  m 


(28) 


The  buoyancy  term  is  equation  (3)  is 


B  -  9.8  2 

6..  P 


(30) 


where 
given  by 


3y  the  local  virtual  potential  temperature  is 


6y  =  (l+.61m)9 


e  is  the  ambient  virtual  temperature. 
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The  water  content  is  given  by 


K 


£ 


is 


m. 


Q  *  ??xfdxdy 
oo  v 


The  condensation  term  in  equation  (5)  is 


G  =  7?Xf . dXdY 
oo  t 


(32) 


(33) 


where  X  is  the  mass  of  the  water,  Y  is  the  mass  of  the  nucleus, 
and  f  is  the  particle  density. 

The  stochastic  coalescence  is  given  by  two  integral  expressions: 

hb x*y’z)  =  l'I!f{Xc*YC*x»y*z)V(Xc’Yc/X,’Y,*x»y*z) 

f(X\Y‘  ,x,y,z)dX‘dY'  (34) 

I(J§|.x,y,z)  *  !!m,Y,x,y,z)V(X,Y/X\Y',x,y,2) 

f (X* ,YJ ,x9y,z)dX'dY‘  (35) 

where  Xc  -  X-X',  Yc  -  Y-Y',  Y  i*>  the  collection  kernel,  and 
the  primed  variables  are  the  parameters  of  the  particle  being 
captured. 


2.1.1  Two  Dimensional  Eulerian  Warm  Tog  Seeding  Model 


The  system  of  equations  governing  this  model  are 


Mt  *  -  ?  ♦  ‘"xx  +  ‘"zz 


ft  •  <V>z  *  K<fxx  +  fzz>  +  m  f>X 


where 


(2506.612M*714)/e 


and  9  is  given. 

These  equations  are  reduced  from  the  general  equations 
given  in  section  2.1,  and  are  valid  only  under  the  follow¬ 
ing  assumptions: 


a>  “t  *  0 

t>)  =  0,  temperature  variations  due  to  condensations 

processes  are  considered  small. 

c>  yw 

d)  and  since  e  is  constant  ever  the  entire  domain, 

than  3K  _  3K  _  n 
then  JX  -  II  '  0 

e)  no  coalescence 


A  schematic  of  the  space  domain  is  represented  in  the  figure 
below  „ 


-7- 


The  boundary  conditions  are: 


a)  No  diffusion  across  the  top,  left  and  bottom  boundaries 
and  symmetry  about  the  vertical  axis. 


ff(x,zL) 

'  I?'*-0) 

0 

(4) 

H<>-v 

*  {£*■•> 

*  go.*)  * 

0 

(5) 

b) 

The  boundary  at  the 

right  is  sui 

Ffi  c 

iently 

far  enough 

away  from 

the  seedii 

ng  area  that 

the 

boundary  conditions 

are: 

M(xl,z) 

*  M,{z) 

(6) 

f(xL,z) 

■  fI 

(7) 

where  Mj  and 

fj  are  the 

initial  com 

tiiti 

ons . 

The 

initial  conditions 

over  the  entire 

domain 

are: 

M(x,z)  * 

Mj  (Z) 

(8} 

-♦» 

X 

* 

NJ 

<S~M* 

II 

fI 

(9) 
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2.1.2  One  Dimensional  Harm  Cumulus  Cloud  Model 


The  reduced  set  of  equations  for  this  model  are  given  as: 


wt  -  -ww  +  Kw  +  8  +  Ew 

£,  £  2 


9t  =  -"e2  *  +  K»zz  +  Ei6 


Mt  =  ‘WHZ  “  p  KKZZ  +  E6M 


ft  ■  -(w+vt)f2  +  Kf22  -  (jff)x  ♦  Ifla1h-Iloss 

where  E,  the  dynamic  entrainment  term  satisfying  two 
dimensional  continuity  is  given  as: 


E'*  0 


E '  <  0 


where 


.61M 

1+.61M 


-  «(los»)2  -w  fe  +  f  -1 

e 

J  • 


where 


8£  is  the  temperature  of  the  environment 


6M  =  M-M, 


i  Mj  the  ambient  mixing  ratio 


se  =  9-o 


,  e,  the  ambient  temperature 
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The  assumption  that  reduce  the  general  set  equations  to 
those  above  are: 


a)  KH=KH*Kp=K 


The  initial  conditions  for  the  system  are 

M(z)  =  «j(z) 

e(z)  =  e-(z) 
f(z)  =  fj 

The  boundary  conditions  for  this  model  are  those  imposed 
by  equation  (4}  in  section  2.1.3  and  the  bottom  boundary. 


f(0)  =  fj 
M(Q)  =  Mj 


2.1.3  Three  Dimensional  Eulerian  Model  of  Warm  Fog 

The  three  dimensional  Eulerian  model  of  warm  fog  modifi¬ 
cation  by  the  addition  of  heat  sources  at  the  ground  level 
are  Identical  to  the  general  set  in  section  2.1  with  the 
assumption  that  the  effects  of  stochastic  coalescence  are 
negligible,  tie  Igai„  =  •  0) 


The  domain  of  the  equations  is  taken  as  rectangular  and 
therefore,  in  a  three-dimensional  system,  there  are  six 
boundary  conditions  to  consider.  The  figure  below  is  a 
schematic  of  the  domain  2 


The  conditions  are  selected  such  that  the  mass  is  flowing 
into  the  left  and  front  hrundaries,  and  flowing  out  of  the 
right,  back  and  top  boundaries.  The  bottom  boundary  is  con¬ 
sidered  the  ground.  The  boundaries  in  which  mass  is  flowing 
out  have  the  boundary  conditions: 


right  face, 
c  ( x  i » y » ..» t } 


C(xL-6x,y,z,t)2 
C( X^-26x,y ,z ,t) 


(i; 


where  C  is  vector  of 


(2) 
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* 


back  face 

C(x,y,  -fiy,z,t}2 

C(x,y,,z,t)  =  - - - - 

C(x,yL“26y,z,t) 


top  face, 
C(x,y,zL,t) 


2 

C(x,y,zL-6z,t) 
C(x,y,z.  -2oZ,t j 

w 


(3) 


(4) 


The  left,  front  and  bottom  faces  are  maintained  at  the 
initial  conditions  of  the  system.  The  initial  conditions 
are  given  as  some  function  of  z  and  are  chosen  such  that 
the  velocities  are  positive  to  ensure  the  proper  flow  out 
of  the  back,  right  and  top  faces  and  in  the  front  and  left 
faces.  The  velocity  at  the  bottom  faca  is  zero  for  all 
time.  Stated  more  precisely. 


C(0,y,z,t)  =  Cj (z) 

(5) 

C(x,0,z,t)  =  C  j ( z ) 

(61 

C(x,y,0,t)  =  Cj(0) 

(7) 

where  C-(z)  is  the  initial  conditions  for  the  entire  system. 


SBSS&BBass, 
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2.2  Method  of  Calculation 


The  equations  presented  for  the  three  models  are  of  the 
same  general  nature  and  are  integrated  numerically  by  an  ex¬ 
plicit  finite  difference  scheme  to  be  discussed  later.  An 
important  consideration  of  the  method  is  the  different  time 
steps  used  in  the  transport  and  condensation  processes.  The 
reason  for  the  different  time  steps  in  the  procedure  are 
only  economic  relative  to  computer  time.  That  is,  if  ail 
the  equations  were  to  be  integrated  with  a  time  step  neces¬ 
sary  to  integrate  the  condensation  processes,  a  large  amount 
of  computer  time  would  be  necessary.  If  the  condensation 
processes  would  be  allowed  to  go  on  independently  during  the 
time  step  of  the  transport  process,  much  computer  time  could 
be  saved.  The  reason  for  this  is  due  to  the  time  saved  in 
storing  and  fetching  quantities  at  the  adjacent  nodes  in  the 
network.  The  condensation  process  is  independent  of  any 
spatial  considerations  during  one  particular  time  step  and 
therefore,  independent  of  neighboring  nodal  processes.  The 
integration  process  used  in  the  time  domain  is  merely  a 
first  term  Taylor  series  expansion 


f ( t+  6 1 )  *  f ( t )  +  5t{f  (t)} 


(1) 


where  6t  is  chosen  so  as  to  maintain  numerical  stability. 

The  problem  presented  in  equations  (1)  through  (7) 
section  2.1  can  be  divided  into  two  main  processes: 

a)  transport  processes 

b)  condensation  processes 

The  transport  processes  occur  at  a  much  lower  rate  than  do 
the  condensation  processes.  Therefore,  to  integrate  both 
processes  at  the  same  6t  would  be  wasteful.  However,  if  the 
two  processes  could  be  integrated  simultaneously  but  at  dif¬ 
ferent  rates,  our  problem  of  waste  would  be  eliminated. 
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For  convenience,  the  problem  can  be  independently  divided 
into  the  two  processes. 


ft  =  T  +  C  (2) 

If  the  stability  criterion  for  transport  (I)  is  St  and  the 
stability  criterion  for  condensation  { C )  is  At  (At-ct), 
then  tne  equation  (l!  could  be  rewritten  as 

f(t*4t)  *  TAt  +  CSt  (3) 

Equation  (3)  could  then  be  applied  as  many  times  as  it 
takes  to  reach  the  time  t+st.  Oaring  the  time  step  St, 

T  is  only  calculated  once,  whereas  C  is  recalculated  at 
each  step  At.  This  amounts  to  a  considerable  savings 
since  most  of  the  computer  time  is  spent  in  calculating 
T.  i-'nr  example,  let  St  =  3At,  then 

f{t+6t/3)  -  TAt  +  C(f(t) }At  +  f(t)  (4) 

f  ( t+;:st/3)  =  TAt  +  C { f  ( t+A t )  > A t  +  f  ( t + a t )  (5) 

*;t+6l)  *  TA t  +  C{ * ( t+2At) }A t  +  f ( t+Z A t )  (6) 


2 . 3  Finite  Difference  Techniques 

The  equations  for  the  three  models  discussed  ie  earlier 
sections  have  In  common  many  terms  of  the  same  general  form. 
Instead  of  giving  the  complete  set  of  finite  difference 
equations  for  all  models,  the  finite  difference  form  for 
each  general  term  will  be  given. 

The  form  of  time  derivatives  are  approximated  by 

ft  =  (*(t+6t)  -  f< t) }/  At  ( 

Other  terms  of  the  general  form  uv^  will  be  approximated 
by  the  following  finite  difference  scheme  by  Crowley^^ 

uvx  =  l/2ui(v.  +  1-v1 _1 }  /  ax  -  (F+-F‘)  /  ax  (; 


and  terms  of  the  general  form  (uv)x  are: 


(uv)„  =  (F"-F+)  /  ax 


where 


■ 

,  .At  +  ,  v 

_Vi+Vh  v  5ui+ru,>_ 


v  =  (vi  +  l+vi>/2‘  y’  =  C  v i _  ?  +  v i ) / 2 


The  terms  of  the  form  (Kux)x  are  approximated  as 


<KVx 


Ki  +  l*Ki)(  VrM-(Ki*Ki-l)(ur,Ji-l) 
2ax2 


The  equations  representing  particle  Effusion  are  con¬ 
tinuous  over  all  X  (mass  classes)  and  over  all  Y  (water 
classes).  These  equations  are  made  discrete  by  a  logarthmic 
transformation  used  by  Berry(2^  The  domain  of  both  X  and  Y 
are  transformed  as  follows: 


0  e3<^2>'Jo 

and 

Yk  =  Yo  *3<k-1)/ko 


Substituting  equations  (7)  and  (8)  into  equation  (7)  of 
section  2.1  gives 


‘t  ufx 


(w+vt)f^,k  +  (K?f:j*k 


X 


+ 


(9) 


where  3c^»0max  and  k=l»kmax«  Equation  (9)  now  represents 
a  discrete  set  of  simultaneous  difference  equations  that 
span  the  domain  of  water  and  mass  classes  of  interest. 
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PDP-8I 

SOFTWARE 

SYSTEM 


3.1  SCH-8 


Program  SCM-8  makes  possible  ths  selective  recovery 
of  data  from  magnetic  tape  cassettes  using  the  DL-622 
1/4"  magnetic  tape  reader.  Provision  is  made  for  checking 
incoming  data  for  accuracy  of  transmission  and  for  notifying 
the  operator,  by  means  of  teletype  messages,  of  parity 
failures  and  program  status. 

Various  checking  features  of  SCH-8  make  possible  the 
recovery  of  data  over  long  periods  of  time  from  specified 
sources  which  yield  information  subject  to  wide  variation 
due  to  meter  scale  or  range  changing.  Checks  are  also 
performed  on  selected  data  to  screen  out  and  correct 
erroneous  data  values  (straddle)  introduced  by  metering 
devices  during  range  transitions. 

Three  averaging  operations  are  provided  for  special 
handling  of  any  one  or  all  of  the  18  input  channels  and 
provision  is  made  for  generation  of  the  standard  deviation 
of  1  or  all  channels  during  a  program  run. 

An  option  code  (E)  is  provided  which,  when  associated 
with  an  input  channel,  produces  a  teletype  listing  of  that 
channels  values  while  not  contributing  to  any  data  averaging 
operation. 

It  should  be  noted  that  the  range  change  checking  and 
straddle  checking  for  the  specified  channels  commence  with 
the  very  first  scan  values  recovered  at  time  of  crogram 
initialization. 

Any  format,  range  change,  or  straddle  parameter 
modifications  must  be  accomplished  prior  to  program 
ini tial i zation . 
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2  Metro-8  and  Dead  Reckoning 

The  programs  Metro-8  (ME-8),  Dead  Reckoning  (DR-8) 
operate  on  DL-820  meteorological  data  as  generated  by 
program  SCH-8.  ME-8  calculates  true  temperature,  true 
air  speed,  Tacan  position  coordinates  with  altitude  cor¬ 
rection  and  wind  direction  and  velocity  from  Ooppier  drift 
measurements.  Wind  level  averaging  of  the  calculated  wind 
velocity  and  direction  is  also  provided.  DR-8  calculates 
dead  reckoning  coordinates  with  or  without  wind  drift 
effects  in  addition  to  true  temperature,  true  air  speed 
and  Tacan  position.  The  programs  are  written  in  FOCAL 
for  use  with  a  PDP-8I  computer  system  having  8K  storage 
and  a  high  speed  reader  (HSR).  A  discussion  of  the  cal¬ 
culations  mentioned  above  follows: 

A.  True  Temperature  and  Air  Speed 

ME-8  and  DR-8  both  require  calculations  of  true 
temperature  T  and  true  air  speed  V  .  The  formulas 

3 

used  are: 

T(°K)  =  TIND*2/3^6 
1+.1992M2 

where  T^p  is  the  measured  (Rosamont)  temperature  in 
°C  and  M  is  the  Mach  number 

Vkn0ts)  '  + 


where  P  is  the  indicated  pressure  in  millibars. 


B.  Tacan  Position 

ME -8  and  DR-8  both  calculate  Tacan  position  as  follows: 
Let  Dy  be  the  direct  line  Tacan  distance,  6j  the  in¬ 
dicated  Tacan  bearing,  eQ  the  magnetic  deviation,  H  the 

altitude  of  the  airplane,  X„  and  the  coordinates  of 

r  o  o 

the  Tacan  station.  The  component  of  the  distance  along 
the  ground  Dq  of  the  airplane  from  the  Tacan  station 
becomes : 


The  RMS  error  in  curve  fitting  to  fourth  order  is  *50. 6  ft 
for  H  in  the  range  1780  ft  (950  mb)  to  23,560  ft  (400  mb). 


I 


r§ 


An  alternative  Method  2  assumes  that  the  Standard 
Atmosphere  obeys  a  power  law  of  the  form 

H  =  a-bPn  ,  n  =  2/7 

From  two  know  points  (H^.Pj)  and  (h‘2,P0)  one  can  then  generate 
the  equation: 

H9-H,  „  n 

H  *  Hl  +  (.  n'"">(pl  ^  1 
H1  ~l2 

An  attached  calculation,  along  with  the  least  squares  program 
and  calculation,  indicates  the  approximation  involved  in  re¬ 
placing  the  2/7  power  by  9/32  (9/32  is  equivalent  to  first 
raising  a  quantity  to  the  9  power,  followed  by  5  successive 
square  roots  to  avoid  the  use  of  logarithm  and  exponential 
functions  in  FOCAL).  The  results  for  F  in  nautical  miles 
are: 

H  =  .53288+2.29121 (6. 98344-Pn)  n=2/7 

H  •  .53288+2.39629(6. 77455-Pn)  n=9/32 

For  P-j  =  900  mb,  H^=3,240  ft,  Po=500  mb,  H2=18,280  ft  both  forms 

of  the  equation  are  good  to  within  about  IS  relative  to  tabu¬ 
lated  data  in  the  range  900  mb  to  about  400  mb. 

C.  Dead  Reckoning  Calculations 


Consider  two  successive  data  points  i-1  and  i  during 
which  time  interval  At.=t. -t.j  -j  the  averaged  true  air  speed 
is  v.  and  the  averaged  true  heading  is  e..  Then  the  incre¬ 
mental  components  of  the  distance  traveled  during  At.  by 
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daad  reckoning  are: 

AXj  a  v^At^sin  e.  ,  Ay.  =  VjAt.cos  e^ 

At  time  t^,  therefore,  the  dead  reckoning  coordinates  are: 

xi  =  xi-l  +  Axi  •  *1  =  *1-1  +A*i 

and  the  flight  path  may  be  generated  by  successive  integra¬ 
tion  over  small  time  intervals. 

The  effect  of  wind  is  to  shift  by  vector  addition  the 
dead  reckoning  path  through  increments  Ax.w  and  Ay*w  caused 
by  a  wind  vector  of  velocity  W.  and  direction  where: 

Ax..w=  -W.-At.sin  ,  Ay^w=  -W^AtjCOS  4. ^ 

Wind  data  is  given  ir>,  the  for":  of  (W^,^)  for  levels  k=l 
through  k=n  corresponding  to  known  pressures  Pk  (or  altitude). 

The  wind  speed  and  direction  at  the  instantaneous  pressure 

P  is  then  found  by  linear  interpolation  of  the  X  and  Y 

components  between  levels  j  and  j+1  such  that  P.^P.^P.^, 

J  i  J+ 1 » 

namely: 


Wx  =  W j (l-Z)sin  VW.  +  12  sin  $j  +  1 
\  -  Kjd-ZJcos  */Wj  +  1Z  cos  »j+, 
AX j W  s  WxAt.  ,  Ay .W*  WyAtj 

where 

z  -  (p-pj)/n>jt)-pj) 
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D .  Approximations  to  the  Mach  Number  and  Arccosine  Function 


Evaluation  of  the  Mach  number  M  is  required  for  the 
computation  of  the  true  temperature  and  true  a':r  speed. 
In  particular  we  calculate: 


=  5  (l+x)n-l  where  x  =  ^ 


since  in  FOCAL  we  cannot  exponentiate  to  a  fractional  power 
of  (2/7)  we  must  instead  calculate: 


=  5  exp{  (2/ 7)1  n ( 1  +x)}  -1 


which  requires  the  use  of  exponential  and  logarithm  functions 
To  avoid  the  use  of  these  functions  a. id  thus  save  locations 
for  variable  storage,  we  use  the  expansion: 


(Hx)n  =  l+nx+Sin|il  ,2  +  n(,n-))(n-2)„3  t. 


which  for  n-2/7  to  fourth  order  in  x  yields 


m2  _  5x  ri.  ,  .  20  2  95  3-, 

H  -  ^  p4-5x+—  x  -  49 


This  approximation  is  accurate  to  .5«  for  O^x^.5 


The  arctangent  function  is  required  in  evaluating  the 
computed  wind  direction  from  Doppler  drift.  Using  the  law 


of  cosines  to  calculate  the  interior  angle  ebetween  the  wind 

1 


speed  vector  and  the  air  speed  vector  yields  cos  (x)=tJ. 
To  find  e,  FOCAL  uses: 


1 


tan  (  /x) 


Bfim«gfiinr 


?I 


32 

t 

-p 

M 


M 


i?£T&HgmS§l&» 


In  an  effort  to  eliminate  the  FOCAL  representation  of 
tan”^x  we  employ  the  following  polynomial  expansion: 

e  -  '/T-x  (aQ+a^x+agX^+a-jX3)  0<x$l 

a0  =  1.57073 
a1  *  -.21211 
a2  =  .07426 
a3  =  -.01873 

_  c 

which  is  accurate  to  an  absolute  error  of  5x10  (.03°) 

in  the  range  O^x^l . 

The  net  effect  of  eliminating  FEXP,  FLOG  and  FATN 
functions  from  FOCAL  is  to  save  some  255  locations,  or 
about  51  variables. 


E.  Determination  of  Mind  Speed  and  Direction  from  Doppler  Drift 


Let  VA  be  the  vector  velocity  of  the  airplane  through  the 
air  directed  at  a  true  angle  from  North,  let  Vg  be  the 
ground  velocity  vector  (Doppler  speed)  making  an  angle 
(Doppler  frift  angle)  with  respect  to  and  let  Vy  be  the 


resulting  wind  velocity  vector  making  an  angle  0y  with 
respect  to  North  which  causes  the  drift.  Solving  the  wind 
triangle  for  Vy  by  the  law  of  cosines  gives: 


»«*  •  VG2  +  »A*  -  2Vs  C05<eD> 


Further  application  of  the  law  of  cosines  gives  the  interior 
anglee  between  and  Vy,  namely: 


,  ,  VG2+VW2-VA2  VG-VAC0S(6D> 
co  s  { “  =  - ~ 

^VGVW  VW 


Finally  the  wind  direction  Oy  is  calculated  as: 


eW  =  °A  +  eD  +  180  1 


where  the  sign  of  $is  the  same  as  the  sign  of  eD< 


F.  Wind  Level 


ME-8  computes  wind  directions  and  velocity  from  Doppler 
drift  for  each  time  point  according  to  section  E.  Because 
of  the  scatter  in  the  calculations  due  to  data  uncertainties 
and  unsteady  flight  conditions,  a  method  of  averaging  the 
wind  speed  and  direction  over  many  data  records  has  been 
instituted  which  can  be  used  to  form  wind  level  averages  for 
input  into  DR-8.  In  particular  the  components  of  the  cal¬ 
culated  wind  are  written: 


Vx,  =  VW  Sin  •  Vy(  '  VK  C0S  »W 


and  added  to  sums  S  and  S  for  pressure  level  i  accordina 

.  x  i  y  i 

to: 


-?6- 


> 


» 

8= 


l 


i 


s 


1 

z 

k=l 


vy 


ik 


where  k=l,2...  n^  is  the  kth  entry  to  the  sums  for  level  i. 

The  pressure  of  wind  levels  are  P1.=940-80  i  (mb)  (i=l,7),  i.e, 
7  levels  in  80  mb  steps  from  860  mb  to  380  mb,  and  a  set  of 
wind  components  is  assumed  to  belong  to  level  i  if  the 
associated  pressure  P  for  which  the  wind  components  were  cal¬ 
culated  falls  within  a  band  P^*40  mb.  At  the  termination  of 
the  run,  the  wind  averages  are  calculated  as: 


V 


1 

i  x 


ei  =  cos-1  (Vy  /91) 


V 


+ 


G •  Dead  Reckoning  Flight  Curvature  Correction 

In  earlier  versions  of  DR-8  dead  reckoning,  no  account 
of  flight  curvature  was  made  other  than  to  assume  that  the 
path  was  straight  between  time  t-j  and  t 2  with  an  average 
heading  of  ^2  and  a  straight  line  distance  travelled  of 
Ds  =  V2(t2-ti )  where  V2  is  the  true  air  speed  calculated  at 
t2.  Consequently,  dead  reckoning  coordinates  were  cal¬ 
culated  as; 


x2  =  Xj+Ds5inij<2  y2  =  y^+Dscos^2  D$  =  V2At 


I 


This  amounts  to  having  the  airplane  flying  in  straight  line 
segments  after  instantaneous  changes  in  heading  angle. 

A  more  realistic  approach  is  to  observe  the  change  in 
heading  and  assume  that  the  flight  path  was  a  circle 

of  radius  R=D/A^{ rad )  where  D  is  the  arc  length  along  the 

V  w 

circle  between  points  1  and  2.  Dc  is  just  Vet 
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and  V=(V-j+V2)/2  is  the  average  velocity.  Consequently  * 
the  straight  line  distance  or  chord  length  D$  between 
points  1  and  2  can  be  shown  to  be: 

D  =f0  f=s  in  (  a  4>  )/a<p  A4>(rad)-n  '4*  ( deg )  / 1 80 

•J  L 


where  f  is  a  fraction  which  changes  from  1.0  for  ao  =  0  to 
.90  for  A^=  90 ° .  Thus  after  a  90°  turn,  the  straight  line 
distance  is  some  10%  smaller  than  that  calculated  with  the 
old  model.  Finally,  the  dead  reckoning  coordinates  of  point 
2  are  now  calculated  as: 

X2=Xi+  2‘^Dc  +  sin*2 1  > 

*"  r  ~i 

Y2  =  Y1  +  lfDc  Jcos^l  +  c° s ’^2  j 


\ 


\ 

I  Pft-m 


More  sophisticated  models  can  be  devised  other  than 
assuming  circular  flight  paths  between  time  elements  when 
angular  changes  occur,  however,  as  a  practical  matter  for 
short  time  intervals,  the  angula"  changes  A<p  will  be  small 
(f  close  to  unity,  5  close  to  and  the  effects  of  flight 
path  curvature  will  be  small  but  cumulative. 
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3.4  FLIPPER 


FLIPPER  is  a  FOCAL  program  for  converting  the  Pencil 

Follower  program  output  data  in  serial  form  to  a  parallel 

form  for  use  as  input  in  other  programs.  In  versions  1.0 

and  1.1  of  Flipper,  the  first  data  string  appearing  on  the 

Pencil  Follower  output  tape  is  treated  as  the  master  curve 

(x!,  yj,  i=l,N<).  Each  succeeding  string  of  data  (x.n»  y.n. 

ja!,N  .n>l)  is  then  processed  by  linear  interpolation  to 
n 

produce  values  cf  ,y  of  that  string  corresponding  to  each 
x.  value  of  the  master  string,  i.e. 


y."  -  y,"  *  {*}-«")  n>,: 

«  J  yn  r»  t  j 


xj+rxo 


3  =  1, N, 


X  n<V  1  < 

xj  xi  xj+l 


Version  1.2  operates  in  a  similar  fashion  except  that 

the  x.’s  are  prescribed  instead  of  saved  from  the  sister 
*  » 

string.  After  the  Pencil  Follower  tape  is  read  with  the 
HSR,  ar  option  is  available  to  select  fast  output  on  the 
HSP  or  a  typed  listing  on  the  TTY.  Version  1.0  has  a  stor¬ 
age  capacity  of  230  data  points,  while  versions  1.1  and  1.2 
can  handle  over  1000  points  with  the  use  of  a  special  FJiEW 
function  which  stores  variables  in  Field  1  of  8K  FOCAL. 


Fog  Drop  Size  Distribution  Program 


r 


i 


0  C 
4/.  J 


This  program  is  designed  to  operate  on  measurements  of 
the  cross  section  of  water  droplets  on  a  photographic  plate 
using  the  curye  digitizer  apparatus.  The  binary  curve  digi¬ 
tizer  output  tape  is  first  processed  by  the  Pencil  Follower 
program  to  rotate  and  translate  the  data  with  reference  to 
the  curve  digitizer  table  and  produce  a  corrected  FOCAL  com¬ 
patible  data  tape.  This  tape  (an  example  of  t'.e  format  is 
attached)  is  now  processed  by  the  Fog  Drop  Size  Distribution 
Program  which  calculates  a  mean  diameter  for  eacii  set  of  4 
edge  measurements  (x^,  y.;  1  =  1  , 4 )  according  to: 


)]  =  a/  (x2-xt)2  +  (y2-y  1 ) 2 


D2  =l/(x4-xj2+(y4-y3)2 


D  =  J  D.  D„ 
mean  1  2 


and  then  produces  a  scaled  histogram  plot  according  to  the 
number  of  classes  and  step  size  specified  by  the  operator. 
The  program  is  written  in  FOCAL,  requiring  8k  extended  stor- 
age,  high  speed  paper  tape  read,  and  no  extended  functions. 


3.6  Vertical  Velocity  Program  for  F100 


The  Vertical  Velocity  program  calculates  the  vertical 
velocity  of  air  from  data  recorded  on  the  F100  oscillograph 
trace.  It  uses  angle  of  attack,  pitch  angle  and  relative 
height  data  reduced  to  the  input  form  via  the  Pencil  Follower 
and  Flipper  (vers  1.2)  programs.  A  mean  True  Air  speed  is 
entered  via  the  teletypewriter  for  each  pass. 

Output  is  via  punched  paper  tape. 

Format  of  output  tape  is: 

Time  (sec) 

Vertical  Velocity,  w,  (m/sec) 

Angle  of  attack,  a,  (deg) 

Pitch  angle,  e*  (deg)  -  with  respect  to  e  at  start 

of  cloud  penetration 

Relative  height,  H,  ( m ) ,  *  with  respect  to  a  value 

prior  to  cloud  penetration. 
Output  is  midpoint  value  of 
3-pcint  averaging  interval. 

Height,  change,  aH,  (m),  -  &H  =  -  H-j 

Angle  of  attack  minus  pitch  angle,  (a  -  <))  (deg) 
Mathematical  description 


1.  (Vertical  velocity  of  air)  -  (vertical  velocity  of  aircraft) 


and  W 


„  aH 

acft  ~  At 


s  Wair  'WacU 


2.  m  sin*  -  TX5  sin  a  *Wair  -  $ 
or  (since  sin«=a  s  1  rve  =e ) 


U-e)  TAS  *  Wa1r  -  AH/at 


I 


Rewriting  3. 


AH 


4.  Uair  *(«-»)  TAS  +  it 


Due  to  noise,  data  are  smoothed  over  a  3-point  array. 
The  smoothing  equations  used  are: 


-  1.  + 


ah 


i  (Hi+rHi-i> 


At 


(T.+,-T.),  (since  time  increments  are  equal.) 

Due  to  possible  precession  of  the  Pitch  Gyro,  a  correction 
term  is  added  to  the  pitch  angle.  This  term  when  applied  to 
initial  data,  (prior  to  cloud  penetration,  i.e.  when  "’'At  is 
small)  forces  the  vertical  velocity  to  be  initially,  approxi¬ 
mately,  zero.  The  term,  TZ  =  (u-T)  is  calculated  using  the 
mean  of  the  first  two  data  sets  available.  To  partially  off¬ 
set  aH  changes,  HT  (1)  is  set  equal  to  HT  (2)  for  initial  data 
sets.  This  is  an  additional  improvement  for  the  calculation 
setting  =  0.  See  Section  7.0  of  program. 


The  program  v.rsion  of  Equation  4  is  then: 
W  =  (AB  -  CTB  +  TZJ)  *TS  *DR  *CT  +  DH/DT 


where  AB 
TB 
TS 
DH 


-  a 


8 

TAS 

AH 


DR 

CT 


IT 

W5 


6080/(3600x3.28)  if  TAS  is  given  in  knots 
1  if  TAS  is  given  in  m/sec 
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3.7  Additional  Software 


The  following  additional  programs  were  developed  to 
provide  a  means  of  operating  and  evaluating  the  performance 
of  various  peripheral  devices  added  to  the  PDP-8I. 

Included  are  specialized  diagnostic  programs,  hardware 
simulation  programs,  and  various  data  handling,  listing,  edit, 
and  code  conversion  programs. 

1 .  CSC 

Column  Selector  Program  provides  capability  of  selecting 
up  to  10  columns  of  data  from  an  input  tape  of  any  consistent 
format  and  reformatting  them  to  be  output  on  either  the  TTY 
or  HSP.  (operator  selection) 

2.  LI BM 

List  IBM  -  Provides  a  means  of  quickly  listing  the  out¬ 
put  of  the  Pencil  Follower  which  is  in  IBM  code  on  the  PDP-8 
Teletype.  (ASCII).  The  output  format  is  defined  by  the  operator 
Pencil  Follower  function  codes  are  interpreted  as  TTY  carriage 
control s . 

3.  fill 

Dump  program  for  the  DL-622.  Reads  buffer  full  of  data 
from  magnetic  tape  and  prints  contents  on  PDP-8  TTY. 

4*  thsr 

Test  High  Speed  reader.  Tests  for  known  sequence  of 
character  and  sends  alarm  if  not  encountered. 

5.  FINTP 

Routine  to  load  core  buffer  from  punched  paper  tape 
containing  specially  formatted  data. 
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6.  DUMP 

Program  to  dump  special  binary  packed  data  onto  paper 

tape . 

7.  RCON 

Program  to  unpack  special  binary  packed  data  from  paper 
tape  and  load  into  field  1  of  PDP-8. 

8.  SIMCYP 

Simulate  operation  of  data  transfer  from  DL-622  to  CIPHER 
utilizing  PCI  and  parity  characters  as  counting  triggers. 
Tranfers  800  characters  of  DL-622  data. 

9.  DGNOS 

IOT  exerciser  fo-  DL-622.  Provides  capability  for 
specifying  delay  in  microseconds  following  "Read  Initiation" 
IOT. 

10.  EDITS 

Tape  copy  with  edit  program.  Copies  punched  paper  tape 
with  character  string  substitution  when  operator  specified 
character  sequence  is  encountered. 

11.  SIMCYP  (Modified) 

DL-622  data  transfer  to  CIPHER  simulator.  Same  as  SIMCYP 
with  the  following  exception:  Counting  triggered  by  occurance 
of  PCI  and  oarity  characters.  A  count  is  maintained  until 
*10'  scans  have  been  detected,  at  which  time  a  'data  block' 
is  transferred. 

12.  SSJ_ 

DL-622  recovery  test  program.  Utilizes  PCI  and  parity 
for  detection  of  BAD  SCANS.  Function  is  to  search  and  log 
bad  scans.  Longitudinal  parity  checking  is  basis  for  identi¬ 
fication. 
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